/*!	 curve_helper.cpp
**	 Curve Helper File
**
**	Copyright (c) 2002-2005 Robert B. Quattlebaum Jr., Adrian Bentley
**
**	This package is free software; you can redistribute it and/or
**	modify it under the terms of the GNU General Public License as
**	published by the Free Software Foundation; either version 2 of
**	the License, or (at your option) any later version.
**
**	This package is distributed in the hope that it will be useful,
**	but WITHOUT ANY WARRANTY; without even the implied warranty of
**	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
**	General Public License for more details.
**
*/

#ifdef USING_PCH
#	include "pch.h"
#else
#ifdef HAVE_CONFIG_H
#	include <config.h>
#endif

#include "curve_helper.h"

#include <algorithm>
#include <vector>

#endif

using namespace std;
using namespace etl;
using namespace synfig;

#define ERR	1e-11
const Real ERROR = 1e-11;

Real synfig::find_closest(const etl::bezier<Point> &curve, const Point &point,
                          float step, Real *dout, float *tout)
{
    Real d, closest = 1.0e50;
    float t, time, closestt = -1;
    Vector p0, p1, end;

    if (dout && *dout > 0) {
        closest = *dout;
    }

    p0 = curve[0];
    end = curve[3];

    for (t = step; t < 1; t += step, p0 = p1) {
        p1 = curve(t);
        d = line_point_distsq(p0, p1, point, time);

        if (d < closest) {
            closest = d;
            closestt = t - step + time * step;
        }
    }

    d = line_point_distsq(p0, end, point, time);

    if (d < closest) {
        closest = d;
        closestt = t - step + time * (1 - t + step); // time between [t-step,1.0]
    }

    // set the time value if we found a closer point
    if (closestt >= 0) {
        if (tout) {
            *tout = closestt;
        }
    }

    return closest;
}

// Line and BezHull Definitions
void BezHull::Bound(const etl::bezier<Point> &b)
{

    // with a starting vertex, find the only vertex that has all other vertices on its right
    int i, j;
    int first, cur, last;

    float d, ds;

    Vector n, vi;
    Vector::value_type	deqn;

    // get left most vertex
    d = b[0][0];
    first = 0;

    for (i = 1; i < 4; ++i) {
        if (b[i][0] < d) {
            d = b[i][0];
            first = i;
        }
    }

    cur = last = first;
    size = 0;

    // find the farthest point with all points on right
    ds = 0;

    do { // should reassign cur so it won't break on first step
        for (i = 0; i < 4; ++i) {
            if (i == cur || i == last) {
                continue;
            }

            // rotate vector to right to make normal
            vi = -(b[i] - b[cur]).perp();
            d = vi.mag_squared();

            // we want only the farthest (solves the case with many points on a line)
            if (d > ds) {
                ds = d;
                deqn = n * b[cur];

                for (j = 0; j < 4; ++j) {
                    d = n * b[i] - deqn;

                    if (d < 0) {
                        break;    // we're on left, nope!
                    }
                }

                // everyone is on right... yay! :)
                if (d >= 0) {
                    // advance point and add last one into hull
                    p[size++] = p[last];
                    last = cur;
                    cur = i;
                }
            }
        }
    } while (cur != first);

}

// Line Intersection
int
synfig::intersect(const Point &p1, const Vector &v1, float &t1,
                  const Point &p2, const Vector &v2, float &t2)
{
    /* Parametric intersection:
    	l1 = p1 + tv1, l2 = p2 + sv2

    	0 = p1+tv1-(p2+sv2)
    	group parameters: sv2 - tv1 = p1-p2

    	^ = transpose
    	invert matrix (on condition det != 0):
    	A[t s]^ = [p1-p2]^

    	A = [-v1 v2]

    	det = v1y.v2x - v1x.v2y

    	if non 0 then A^-1 = invdet * | v2y -v2x |
    								  | v1y -v1x |

    	[t s]^ = A^-1 [p1-p2]^
    */

    Vector::value_type det = v1[1] * v2[0] - v1[0] * v2[1];

    // is determinant valid?
    if (det > ERR || det < -ERR) {
        Vector p_p = p1 - p2;

        det = 1 / det;

        t1 = det * (v2[1] * p_p[0] - v2[0] * p_p[1]);
        t2 = det * (v1[1] * p_p[0] - v1[0] * p_p[1]);

        return 1;
    }

    return 0;
}

// Returns the true or false intersection of a rectangle and a line
int intersect(const Rect &r, const Point &p, const Vector &v)
{
    float t[4] = {0};

    /*get horizontal intersections and then vertical intersections
    	and intersect them

    	Vertical planes - n = (1,0)
    	Horizontal planes - n = (0,1)

    	so if we are solving for ray with implicit line
    */

    // solve horizontal
    if (v[0] > ERR || v[0] < -ERR) {
        // solve for t0, t1
        t[0] = (r.minx - p[0]) / v[0];
        t[1] = (r.maxx - p[0]) / v[0];
    } else {
        return (int)(p[1] >= r.miny && p[1] <= r.maxy);
    }

    // solve vertical
    if (v[1] > ERR || v[1] < -ERR) {
        // solve for t0, t1
        t[2] = (r.miny - p[1]) / v[1];
        t[3] = (r.maxy - p[1]) / v[1];
    } else {
        return (int)(p[0] >= r.minx && p[0] <= r.maxx);
    }

    return (int)(t[0] <= t[3] && t[1] >= t[2]);
}

int synfig::intersect(const Rect &r, const Point &p)
{
    return (p[1] < r.maxy && p[1] > r.miny) && p[0] > r.minx;
}

// returns 0 or 1 for true or false number of intersections of a ray with a bezier convex hull
int intersect(const BezHull &bh, const Point &p, const Vector &v)
{
    float mint = 0, maxt = 1e20;

    // polygon clipping
    Vector n;
    Vector::value_type	nv;

    Point last = bh.p[3];

    for (int i = 0; i < bh.size; ++i) {
        n = (bh.p[i] - last).perp(); // rotate 90 deg.

        /*
        	since rotated left
        	if n.v 	< 0 - going in
        			> 0 - going out
        			= 0 - parallel
        */
        nv = n * v;

        // going OUT
        if (nv > ERR) {
            maxt = min(maxt, (float)((n * (p - last)) / nv));
        } else if (nv < -ERR) { // going IN
            mint = max(mint, (float)((n * (p - last)) / nv));
        } else {
            if (n * (p - last) > 0) { // outside entirely
                return 0;
            }
        }

        last = bh.p[i];
    }

    return 0;
}

int Clip(const Rect &r, const Point &p1, const Point &p2, Point *op1, Point *op2)
{
    float t1 = 0, t2 = 1;
    Vector v = p2 - p1;

    /*get horizontal intersections and then vertical intersections
    	and intersect them

    	Vertical planes - n = (1,0)
    	Horizontal planes - n = (0,1)

    	so if we are solving for ray with implicit line
    */

    // solve horizontal
    if (v[0] > ERR || v[0] < -ERR) {
        // solve for t0, t1
        float 	tt1 = (r.minx - p1[0]) / v[0],
                tt2 = (r.maxx - p1[0]) / v[0];

        // line in positive direction (normal comparisons
        if (tt1 < tt2) {
            t1 = max(t1, tt1);
            t2 = min(t2, tt2);
        } else {
            t1 = max(t1, tt2);
            t2 = min(t2, tt1);
        }
    } else {
        if (p1[1] < r.miny || p1[1] > r.maxy) {
            return 0;
        }
    }

    // solve vertical
    if (v[1] > ERR || v[1] < -ERR) {
        // solve for t0, t1
        float 	tt1 = (r.miny - p1[1]) / v[1],
                tt2 = (r.maxy - p1[1]) / v[1];

        // line in positive direction (normal comparisons
        if (tt1 < tt2) {
            t1 = max(t1, tt1);
            t2 = min(t2, tt2);
        } else {
            t1 = max(t1, tt2);
            t2 = min(t2, tt1);
        }
    } else {
        if (p1[0] < r.minx || p1[0] > r.maxx) {
            return 0;
        }
    }

    if (op1) {
        *op1 = p1 + v * t1;
    }

    if (op2) {
        *op2 = p1 + v * t2;
    }

    return 1;
}

static void clean_bez(const bezier<Point> &b, bezier<Point> &out)
{
    bezier<Point> temp;

    temp = b;
    temp.set_r(0);
    temp.set_s(1);

    if (b.get_r() != 0) {
        temp.subdivide(0, &temp, b.get_r());
    }

    if (b.get_s() != 1) {
        temp.subdivide(&temp, 0, b.get_s());
    }

    out = temp;
}

// CIntersect Definitions

CIntersect::CIntersect()
    : max_depth(10) // depth of 10 means timevalue parameters will have an approx. error bound of 2^-10
{
}

struct CIntersect::SCurve {
    bezier<Point> 	b;		// the current subdivided curve
    float rt, st;
    // float 			mid,	// the midpoint time value on this section of the subdivided curve
    //				scale;	// the current delta in time values this curve would be on original curve

    float 	mag;			// approximate sum of magnitudes of each edge of control polygon
    Rect	aabb;			// Axis Aligned Bounding Box for quick (albeit less accurate) collision

    SCurve(): b(), rt(), st(), mag() {}

    SCurve(const bezier<Point> &c, float rin, float sin)
        : b(c), rt(rin), st(sin), mag(1)
    {
        Bound(aabb, b);
    }

    void Split(SCurve &l, SCurve &r) const
    {
        b.subdivide(&l.b, &r.b);

        l.rt = rt;
        r.st = st;
        l.st = r.rt = (rt + st) / 2;

        Bound(l.aabb, l.b);
        Bound(r.aabb, r.b);
    }
};

// Curve to the left of point test
static int recurse_intersect(const CIntersect::SCurve &b, const Point &p1, int depthleft = 10)
{
    // reject when the line does not intersect the bounding box
    if (!intersect(b.aabb, p1)) {
        return 0;
    }

    // accept curves (and perform super detailed check for intersections)
    // if the values are below tolerance

    // NOTE FOR BETTERING OF ALGORITHM: SHOULD ALSO/IN-PLACE-OF CHECK MAGNITUDE OF EDGES (or approximate)
    if (depthleft <= 0) {
        // NOTE FOR IMPROVEMENT: Polish roots based on original curve
        //						(may be too expensive to be effective)
        int turn = 0;

        for (int i = 0; i < 3; ++i) {
            // intersect line segments

            // solve for the y_value
            Vector v = b.b[i + 1] - b.b[i];

            if (v[1] > ERROR && v[1] < ERROR) {
                Real xi = (p1[1] - b.b[i][1]) / v[1];

                // and add in the turn (up or down) if it's valid
                if (xi < p1[0]) {
                    turn += (v[1] > 0) ? 1 : -1;
                }
            }
        }

        return turn;
    }

    // subdivide the curve and continue
    CIntersect::SCurve l1, r1;
    b.Split(l1, r1);	// subdivide left

    // test each subdivision against the point
    return recurse_intersect(l1, p1) + recurse_intersect(r1, p1);
}

int intersect(const bezier<Point> &b, const Point &p)
{
    CIntersect::SCurve	sb;
    clean_bez(b, sb.b);

    sb.rt = 0;
    sb.st = 1;
    sb.mag = 1;
    Bound(sb.aabb, sb.b);

    return recurse_intersect(sb, p);
}

// Curve curve intersection
void CIntersect::recurse_intersect(const SCurve &left, const SCurve &right, int depth)
{
    // reject curves that do not overlap with bounding boxes
    if (!intersect(left.aabb, right.aabb)) {
        return;
    }

    // accept curves (and perform super detailed check for intersections)
    // if the values are below tolerance

    // NOTE FOR BETTERING OF ALGORITHM: SHOULD ALSO/IN-PLACE-OF CHECK MAGNITUDE OF EDGES (or approximate)
    if (depth >= max_depth) {
        // NOTE FOR IMPROVEMENT: Polish roots based on original curve with the Jacobian
        //						(may be too expensive to be effective)

        // perform root approximation
        // collide line segments

        float t, s;

        for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < 3; ++j) {
                // intersect line segments
                if (intersect_line_segments(left.b[i], left.b[i + 1], t, right.b[j], right.b[j + 1], s)) {
                    // We got one Jimmy
                    times.push_back(intersect_set::value_type(t, s));
                }
            }
        }

        return;
    }

    // NOTE FOR IMPROVEMENT: only subdivide one curve and choose the one that has
    //						the highest approximated length
    // fast approximation to curve length may be hard (accurate would
    // involve 3 square roots), could sum the squares which would be
    // quick but inaccurate

    SCurve l1, r1, l2, r2;
    left.Split(l1, r1);	// subdivide left
    right.Split(l2, r2); // subdivide right

    // Test each candidate against each other
    recurse_intersect(l1, l2);
    recurse_intersect(l1, r2);
    recurse_intersect(r1, l2);
    recurse_intersect(r1, r2);
}

bool CIntersect::operator()(const etl::bezier<Point> &c1, const etl::bezier<Point> &c2)
{
    times.clear();

    // need to subdivide and check recursive bounding regions against each other
    // so track a list of dirty curves and compare compare compare

    // temporary curves for subdivision
    CIntersect			intersector;
    CIntersect::SCurve	left, right;

    // Make sure the parameters are normalized (so we don't compare unwanted parts of the curves,
    //	and don't miss any for that matter)

    // left curve
    // Compile information about curve
    clean_bez(c1, left.b);
    left.rt = 0;
    left.st = 1;
    Bound(left.aabb, left.b);

    // right curve
    // Compile information about right curve
    clean_bez(c2, right.b);
    right.rt = 0;
    right.st = 1;
    Bound(right.aabb, right.b);

    // Perform Curve intersection
    intersector.recurse_intersect(left, right);

    // Get information about roots (yay! :P)
    return times.size() != 0;
}

// point inside curve - return +/- hit up or down edge
int intersect_scurve(const CIntersect::SCurve &b, const Point &p)
{
    // initial reject/approve etc.


    if (p[0] < b.aabb.minx || p[1] < b.aabb.miny || p[1] > b.aabb.maxy) {
        return 0;
    }

    // approve only if to the right of rect around 2 end points
    {
        Rect 	r;
        r.set_point(b.b[0][0], b.b[0][1]);
        r.expand(b.b[3][0], b.b[3][1]);

        if (p[0] >= r.maxx && p[1] <= r.maxy && p[1] >= r.miny) {
            float df = b.b[3][1] - b.b[0][1];

            return df >= 0 ? 1 : -1;
        }
    }

    // subdivide and check again!
    CIntersect::SCurve	l, r;
    b.Split(l, r);
    return	intersect_scurve(l, p) + intersect_scurve(r, p);
}

int synfig::intersect(const bezier<Point> &b, const Point &p)
{
    CIntersect::SCurve	c(b, 0, 1);

    return intersect_scurve(c, p);
}